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ABSTRACT 



In this paper, we further develop the model for the migration of planets intro- 
duced in Del Popolo et al. (2001). We first model the protoplanetary nebula as a 
time-dependent accretion disc and find self-similar solutions to the equations of the 
accretion disc that give to us explicit formulas for the spatial structure and the tempo- 
ral evolution of the nebula. These equations are then used to obtain the migration rate 
of the planet in the planetesimal disc and to study how the migration rate depends on 
the disc mass, on its time evolution and on some values of the dimensionless viscosity 
parameter a. We find that planets that are embedded in planetesimal discs, having 
total mass of 10^'* — O.lAf0, can migrate inward a large distance for low values of 
a (e.g., a ~ 10^'^ — 10^^) and/or large disc mass and can survive only if the inner 
disc is truncated or because of tidal interaction with the star. Orbits with larger a are 
obtained for smaller value of the disc mass and/or for larger values of a. This model 
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may explain several orbital features of the recently discovered giant planets orbiting 
nearby stars. 

Key words: Planets and satellites: general; planetary system 

1 INTRODUCTION 

After the discovery of 51 Peg by Mayor & Queloz (1995), more than sixty extrasolar planet candidates have been discovered. 
These planets have unexpected properties: three planets (51 Peg, r Boo, v And) are in extremely tight circular orbits with 
periods of a few days, two planets (p^ Cue and p CrB) have circular orbits with periods of order tens of days and three planets 
with wider orbits (16 Cyg B, 70 Vir and HD 114762) have very large eccentricities. Between the unexpected properties of 
these planets, most of which are Jupiter-mass objects, particularly noteworthy is the small orbital separations at which these 
planets orbit around their parent stars: among the several tens of planets detected so far, at least fifteen of the planets orbit 
at a distance between ~ 0.046 and O.llAU from their parent star. The properties of these planets are difficult to explain using 
the quoted standard model for planet formation (Lissauer 1993; Boss 1995). This standard model predicts nearly circular 
planetary orbits, and giant orbital distances > 1 AU from the central star so that the temperature in the protostellar nebula 
is low enough for icy materials to condense (Boss 1995, 1996; Wuchterl 1993, 1996). 

The most natural explanation for planets on very short orbits is that these planets have formed further away in the 
protoplanetary nebula and they have migrated to the small orbital distances at which they are observed. Some authors have 
also proposed scenarios in which migration and formation were concurrent (Terquem et al. 1999). 

So far, four mechanisms have been proposed to explain the presence of planets at small orbital distances. A first mechanism 
deals with dynamical instabilities in a system of giant planets (Rasio & Ford 1996; Weidenschilling & Marzari 1996). The orbits 
of planets could become unstable if the orbital radii evolve secularly at different rates or if the masses increase significantly as 
the planets accrete their gaseous envelopes (Lissauer 1993). In this model, the gravitational interaction between two planets, 
during evolution, (Gladman 1993; Chambers et al. 1996) can give rise to the ejection of one planet, leaving the other in 
a smaller orbit. While it is almost certain that this mechanism operates in many systems with multiple planets, it cannot 
account for the relatively large number of short-period planets observed (Terquem et al. 1999). 

The second mechanism, called 'migration instability' (Murray et al. 1998), involves a resonant interaction between the planet 
and a disc of planetesimals, located in its orbit which leads to the planetesimals ejection and the inward migration of the 
planet. The advantage of this mechanism is that the migration is halted naturally at short distances when the majority of 
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perturbed planetesimals collide with the star. Moreover wide eccentric orbits can also be produced for planets more massive 
than ~ 3Afj. However the model has some disadvantages, since the protoplanetary disc mass required for the migration of a 
Jupiter-mass planet to a ~ 0.1 AU is very large (Ford et al. 1999; Terquem et al. 1999). ^ 

The third possible mechanism proposed to explain short period of the planets is the dissipation in the protostellar nebula 
(Goldreich & Tremaine 1979, 1980; Ward 1986; Lin et al. 1996; Ward 1997). Since, in this model, the time-scale of migration 
is ~ lO^xp-yr (Ward 1997), the migration has to switch off at a critical moment, if the planet has to stop close to the star 
without falling in it. The movement of the planet might be halted by short-range tidal or magnetic effects from the central 
star (Lin et al. 1996) (in any case, as shown by Murray et al. (1998), it is difficult to explain, by means of these stopping 
mechanisms, planets with semi- major axes a > 0.2 AU). 

The fourth mechanism is based upon dynamical friction between the planet and a planetesimal disc (Del Popolo et 
al. 2001) ^. In that paper, we showed that dynamical friction between a planet and a planetesimals disc is an important 
mechanism for planet migration and showed that migration of a IMj planet to small heliocentric distances (0.05 AU) is 
possible for a disc with a total mass of lO"'' -i- W^^Mq if the planetesimal disc does not dissipate during the planet migration 
or if the disc has Md > O.OIMq and the planetesimals are dissipated in ^ 10* yr. The model predicts that massive planets 
can be present at any heliocentric distances for the right value of disc mass and time evolution. 
Some advantages of the model are: 

1) differently from models based on the density wave theory (Goldreich & Tremaine 1980; Ward 1986, 1997), our model: 

a) does not require a peculiar mechanism to stop the inward migration (Lin et al. 1996). Planet halt is naturally provided by 
the model. 

b) It can explain planets found at heliocentric distances of > 0.1 AU or planets having larger values of eccentricity. 

c) It can explain metallicity enhancements observed in stars having planets in short-period orbits. 

2) Whereas the model of Murray et al. (1998) has the drawback of requiring very massive discs (Ford et al. 1999; Terquem et 
al. 1999), our model shows that radial migration is possible with modest masses of planetesimals discs and predicts the right 
metallicity enhancement. 

A point that requires improvement in our model is the model used for the planetesimal disc. In fact, in Del Popolo et al. (2001), 
we used Opik (1976) approximation, assuming that the surface density in planetesimals varies as Ss(»') = ^©(lAU/r)'''''^, 

* The migration of a Jupiter mass planet from 5 AU to very small radii requires about O.IMq in 5 AU. 

t In contrast to the hypothesis proposed by Fernandez &; Ip (1984), migration in the model in Del Popolo et al. (2001) does not require 
the presence of asymmetry in the planetesimals distribution 
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where E©, the surface density at 1 AU, is a free parameter. The model described in Del Popolo et al.(2001) can be improved 
using a more reliable model for the disc, and in particular using a time-dependent accretion disc, since it is widely accepted 
that the solar system at early phases in its evolution, is well described by this kind of structure. Moreover, astronomical 
observations of the last decade have led to the conclusion that discs around young stellar objects, for example T Tauri stars, 
are Keplerian accretion discs and that they are ever-changing and having a limited life-span. In the next section, we shall get 
self-similar solutions of the diffusion equation and will build up analytical explicit formulas for the surface density and the 
other physical quantities required to calculate the planet migration, due to dynamical friction between the planet and the 
planetesimals in the disc. 

The plan of the paper is the following: in Sect. 2, we introduce the disc model used to study radial migration. In Sect. 
3, we review the migration model introduced in Del Popolo et al. (2001). In Sect. 4 we show the results that can be drawn 
from our calculations and finally the Sect. 5 is devoted to the conclusions. 

2 DISC MODEL 

In order to study the formation of planetary systems it is necessary to study the global evolution of solid material which 
constitutes, together with gas, the protoplanetary discs. The idea that discs have an important role in stars and planet 
formation is not a new one: papers by Peek (1942), von Weizsaker (1943,1948) and Liist (1952) introduced the idea that 
the solar nebula was an accretion disc, while the seminal paper by Lynden-Bell & Pringle (1974) foreshadowed the present 
view that discs are commonly found in early stellar formation (Beckwith et al. 1990). In fact, protostellar discs around 
young stellar objects that have properties similar to that supposed for the solar nebula are common: between 25 to 75% of 
young stellar objects in the Orion nebula seem to have discs (Prosser et al. 1994; McCaughrean & Stauffer 1994) with mass 
1O"^M0 < Md < 1O"^M0 and size 40 ±20 AU (Beckwith & Sargent 1996). The previous quoted evidences had led to a large 
consensus about the nebular origin of the Solar System. Moreover, the observations of circumstellar discs surrounding T Tauri 
stars support the view of a disc having a limited life-span and characterized by continuous changes during its life. This means 
that models like the minimum-mass model cannot model properly the solar nebula: in order to model the spatial and temporal 
changes of the disc, one must use the theory of time-dependent accretion disc. While numerical models of time-dependent 
accretion disc (see Ruden & Lin 1986; Ruden & Pollack 1991; Reyes-Ruiz & Stepinski 1995) are used in the astronomical 
context, in the cosmogonic one they have not been widely used, at least till recent. In the past, but also in several recent works, 
the origin of the Solar System was based upon a steady-state model, namely the minimum-mass model (Weidenschilling 1977; 
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Hayashi 1985). Other papers assumed only spatial, but not temporal, changes in the nebula and the origin and evolution of 
the Solar System was studied by means of steady-state accretion disc models (see MorfiU & Wood 1989; Stepinski et al. 1993). 
For the reasons previously described and for others stressed by Stepinski (1998), the temporal evolution of the nebula must 
be taken into account, and this can be accomplished by means of the time-dependent accretion disc model. In the following, 
we introduce a time- dependent accretion disc model that shall be used in the next sections to study planets migration. 
Before starting it is useful to divide, as customary, the evolution of the solar nebula in three stages: 

a) the formation stage, in which the nebula is build up by infalling matter; 

b) the viscous stage, in which internal torques produce the redistribution of angular momentum; 

c) the clearing stage, in which the gaseous component of the nebula is dispersed. 

The equation that are to be solved to calculate the properties of the nebula are the thin-disc set of equations (Frank et 
al. 1985; Stepinski 1998): 

E = 2Hp (1) 
H = (2) 



2 k 



c, = 



■T (3) 



16ct 4 9 2 

V — acsH (5) 
K = nap^T^ (6) 



as _ 3_a 

dt r dr 



1/2 d . 1/2, 

or 



(7) 



where p is the volume density, E is the surface density, H is the half thickness of the disc, Cs is the sound velocity, 57k = 
GM, I r^ is the Keplerian angular velocity, T is the temperature at the vertical center of the disc, v is the viscosity and k 
is the Rosseland mean opacity, k is the Boltzmann's constant, \i is the mean molecular weight (assumed to be 2.33 (Ruden 
& Pollack 1991)), TTip is the mass of the proton, a is the Stephan-Boltzmann constant and a is the dimensionless viscosity 
introduced by Shakura & Sunyaev (1973). In some papers, it is assumed a = 0.01 as a fiducial value for a disc driven by 
thermal convection (Ruden & Pollack 1991; Reyes-Ruiz & Stepinski 1995). In the following, we shall perform the calculations 
for three values of q: 10~^, 10~^ and 10"^, which span the range from inefficient to highly efficient turbulent convection. See 
that all the radial dependence is carried in Q.k- The previous equations constitute a complete physical model from which 
it is possible to calculate the physical quantities of the nebula. Similar to Stepinski (1998), we have assumed the piecewise 
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power-law formula describing the Rosseland mean opacity, k = n{p,T) given in Ruden & Pollack (1991). The formula we 
used for f — !/(E,r) is that of Reyes-Ruiz and Stepinski (1995). Tcrit represents the critical midplane optical depth required 
for convective viscosity to be present. The value assumed is that of Ruden & Pollack (1991), Tcrit — 1.78 (see also Ruden & 
Pollack (1991) and Reyes-Ruiz & Stepinski (1995) for a discussion concerning the introduction of the critical optical depth). 

In order to calculate the properties of the nebula, we shall take advantage of the fact that the nonlinear diffusion process 
is self-similar. We recall that self-similar solutions in accretion discs have been studied in general by Pringle (1974), Filipov 
(1984), Filipov, Lyubarskii & Shakura (1987); in the context of AGN by Cannizzo, Lee & Goodman (1990); tn the context 
of Dwarf Nova outbursts by Mineshige (1991); m the context of fallback discs around newborn neutron stars by Mineshige, 
Nomoto & Shigeyama (1993) and m the context of Anomalous X-Ray Pulsars by Perna, Hernquist & Narayan (1999). 

Self-similar solutions, asymptotically exact for real problems, give very important information about the nature of vis- 
cously evolving discs, which does not follow directly from the numerical calculations. On the other hand, self-similar disc 
solutions have the unphysical property that the inner edge of the disc (inner radius) is at the origin {Ri„ — 0) , not at the 
surface of the accreting star or somewhere above, like the magnetic radius. 

When the viscous stage starts, namely when T,{r,t) evolution is governed by equation (Q), the nebula has the initial 
distribution left by the formation stage processes. The solution of the nonlinear diffusion equation (equation (Q)) becomes 
self-similar (Ruden 1993) after the short transitional time needed to accommodate the initial conditions and then the self- 
similar evolution can describe the nebula evolution except the transition phase. The latter is only important for the link it 
provides to initial conditions. Since the conditions at the beginning of the viscous stage are poorly constrained, it is convenient, 
as remarked by Stepinski (1998), to assume as initial conditions the distribution of surface density at the beginning of the 
viscous regime. ^ In order to get the solution for E(r, t), we first solve the first five equations algebraically. As an intermediate 
step towards this we solve the equations (|^),(^,(|^,(^) and (^ in terms of Qk, E and T: 




V2a 



k 1 



T 



V 



finip Qk 



t The specific form of initial conditions does not much influence the evolution of the gas, inasmuch as the process is difl^usive in nature 



and the initial conditions are forgotten after a time short in comparison with the evolutionary time scales. 
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p = 



:) 



( — 

\fimp J 

Then we place these results into equation (El) and solve for T in terms of fix and E: 



\ o ) \fimp J 



Now we place this into the equation for u and obtain it in terms of CIk and E: 
, ^ 72^™+^ (27 X 2-^^l ™ E^^™-^ 



Now using Qk ~ \/ GM* / r'-^ we obtain the viscosity in terms of r and E 



^a6-2Va+i 1^27 X 2" 



(GM) '3-2'6+a 2 J- 3(g_2"'^^„ 2)j]l3~2b+' 



which is in the form 



V = Cr^T." 



(8) 



where 



„ /— ^ I 1 / 11 + 3q K.n \ 6-2b+c 



f k \ <5-2i'+f a+1 1 
(GM) + » "2 

\ /imp 



a+1 



6 - 26 + a 2 
2 (a + 2) 



6 - 2fe + a 

The solution of the equation for such a form of the viscosity can be obtained following Pringle (1974), and Mineshige et 
al. (1993), and is given by: 



Eo 
where 

K = 



K 



I r 



^ \Rit)) 



2q-p + 2 ' 
9+1 



(9) 



2q 



(5g-2p + 4)(2g-p + 2)^ 

and the outer radius of the disc, R{t), is: 
, > 2 

t \ 5g-2p+4 



(10) 



The scales ro, to and Eo are free except that they should satisfy the equation: 



(11) 
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Note that for q = the diffusion equation is Unear and the solution is as foUows: 



Md (0) 1 
2-K ilCa^~H 



Q-4i 




( 



) 



S(r,t) 




(12) 



where Md(0) is the initial mass of the disc, I = 2 — p, = 4-2p ^'^'^ -^"^i modified Bessel function of order i^i. Here the 

initial form of E is chosen to be a Dirac-delta distribution at radius a. 



The solution for E(r, t) (equation (|9|)) has an explicit radial dependence while the temporal dependence comes from the 
time dependence of the outer radius. Three scaling parameters are present: ro, Eo, and to, and this last implicitly depends on 
the viscosity. Note that the similarity solutions have four parameters: the initial disc mass, Md, the initial disc characteristic 
radius, Rin, the value of the viscosity and its radial dependence. The solutions must be continuous at the boundaries between 
different viscosity regimes. We indicate with the boundary between the i-th and j-th viscosity regime. In order to calculate 
the boundaries, rij, and impose that E(r, t)i are continuous on them, we may follow the same technique used by Stepinski 
(1998). Namely, we write the surface density in the i-th regime as: 



where Fo,i are constants. At the boundary between the 2 and the 3 regime the temperature is 150 K. We can substitute this 
temperature to the equation for the temperature in the regime 2, T2, to obtain an equation with two unknowns, r23 and the 
factor Fo,2, which comes from the dependency of the temperature on S*. Assuming as normalization Fo,2 = 1, '"23 can be 
calculated solving the equation T2{r,t) — 150K. Once the value of r23 is known, we can calculate Fo,3 in two ways: 

1) using the equation for the temperature for the 3 regime, substituting in it r = r23 and solving T3{r23, Fo^3,t) = 150 for 
Fo,3. 

2) Imposing the condition S2(r23) = E3(r23, -Fo,3) and solving for -Fo,3- 

This procedure can be repeated for regimes 4 and 5. To find the ri2 boundary, we follow the same procedure previously 
used but with optical depth rather than temperature, namely we use the equation for the opacity, t = ^ (Ruden & Pollack 
1991). In this way, it is possible to obtain the boundaries between opacity regimes, in the same way as Stepinski (1998) (see 
their equation (15), (16)). Another way of obtaining approximated values for the boundaries is that of approximating the 
temperature using its dependence in the SIGOR regime and solving the same equations. 



4 



4 




E*(r,t) = Fo,« X Ei(r,t) 



(13) 
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Table 1. Summary of opacity and viscosity regimes 



Number 


Regime 


Applicability 


Opacity 


"mf.' 




c 




P 


1? 


1 


MOTOR 




K = 2 X io-5r2 


6 6 2 

0.14a B/i stc-j.™, 5 




6 

5 





2 


IGOR 


T < 150 K 


K = 2 X 10-5t2 


2.03 X 10i°a2/i-2 







2 


3 


IGSOR 


150 K< T < 180 K 


K = 1.15 X 10"T-8 


2.93 X IO-^qTT /^-TT m~ 


5 

TT 


15 
11 


2 
11 


4 


SIGOR 


180 K< T < 1380 K 


K = 2.13 X 10 




13 13 5 

142. a 9 /X 9 m is 




5 
6 


8 
9 


5 


SIGSOR 


T > 1380 K 


K = 4.38 X 1043T-" 


, IS 18 

6. 44 X 10~-^oT7/i T7m 


S 

T7 


24 
17 


2 

17 



Note: m = M,/Mq 
2.0.1 Different Opacity Regimes 

Different opacity regimes witli different ko, cl and h will give different sets of C, p and g. In Table 1, we have summarized 
the opacity and viscosity regimes. As in Stepinski (1998), we have used the following opacity regimes: MOTOR (Marginally 
optically thick opacity regime) (is a subset of IGOR but is based on optical depth instead of temperature) ; IGOR (Ice grains 
opacity regime); IGSOR (Ice grains sublimate opacity regime); SIGOR (Silicate and iron grains opacity regime); SIGSOR 
(Silicate and iron grains sublimate opacity regime). In the case of MOTOR, since 5 = the diffusion equation is linear and 
finding a solution, which can be expressed in terms of Bessel functions, is much easier (see equation (|l^)). It is interesting 
to note the role of opacity in changing the disc properties. For example, an increase in opacity forces the disc to become 
convective much nearer the surface, or in other terms produces a thickening of the convective zone and an increase in the 
convective velocity. This produces an increase in the viscosity leads to a final decrease in the surface density. 

2.1 Determining The Scale Constants 

Integrating the solution, equation (^, gives the mass of the disc as a function of time 

Md= 2nr-Edr^K''+^{5q~2p + 4)TvrlEo(—] (14) 



Let Td be the time when the all effects due to the initial conditions had vanished and the epoch of self-similar evolution has 
started. Td is of the order of dynamical timescale at the inner radius, Rin ^: 



niiiRin) 

If the mass of the disc at Td is M^, then 



M°d = Md{Td) = K'+'{5q -2p + 4)nrlJ^o (^) 

§ In the following wc clioose _Ri„ = 0.037AU. 
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Table 2. Scale constants in different opacity regimes 



Number Regime tq to 



66— 22i /"^ M° \ 5 

MOTOR 1 2.38a-sij,5T^^'im5r^ ( s^Tk- ) (SC^d)^ 



IGOR 



1 1.64xlO-"a-v(TSr)"'(3CT,)-^r,I ^ (3^,)!^ 

IGSOR 1 lU.a-^,^m^^^y^i3CT,r^r^ (_^)«(3CT,)^ 

, 13 13 5 / 9M? \ ~tW , 8 SI / 9J\/0 \ ft , ,9 

SIGOR 1 2.35 X lO-^a" — ^— mis ^ {3CTd)"^ r^** V (3CTd)?^ 

18 18 A / ITA/S , „^ ,_J_ t4 / 17A/5 \ , ,12 

rf^i-mrf <y ] {3CTa) r^^ ^rr (3CTd)^ 



SIGSOR 1 51. So- 



using equation (|ll| ) in this equation, and solving for Eo, we get: 

Eo=^ro"^ (15) 



where 

5g-2p+4 

/ MS \ 1 

^=[ i,,-2p+\)nK.^^ ) (3C^^^)™ (16) 

Using equation ( [l5| ) in equation (pj|), one obtains 

1 5g+4-2p 

io = ^A-% ^ (17) 



Using this last in equation (|10[), we get 
i?(t)=(^A-)"™i™ 
So i?(t) and thus the solutions are independent of the numerical parameter ro. So for simplicity we choose 



ro = 1 

Then we obtain 
1 

to = 



The previous calculations are summarized in Table 2 for different opacity regimes, namely the table gives information on the 
scale parameters ro, to, and Eq. Using D,k = \J GMt, / in equation (^) temperature distribution can be written as 



2-a 



/ ll + 3a \ e-2b + a / fe \ 0-^" + " , ^ (°+l) -3(a + l) 2(a + 2) 

(^27x2 (GM,)^^^ r 6-2i.+a ss^^m^ (18) 
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Equation (^8|) displays the radial dependence of the solutions, the temporal one in terms of the outer radius and the parameters 
Eo, to, Q and jj,. As in the case of E(r, t), the temperature is continuous at nj. 



3 REVIEW OF THE PLANETS MIGRATION MODEL 

The model used to study the migration of extra-solar planets was introduced in two recent papers (Del Popolo et al. 1999, Del 
Popolo et al. 2001). The same model shall be used in the present paper to study the radial migration of extrasolar planets. 
The difference between the present paper and the previous one, regarding the migration of planets in planetesimal discs (Del 
Popolo et al. 2001), is only due to the model used to describe the spatial structure and temporal evolution of the disc. In the 
following, we review the planets migration model introduced in Del Popolo et al. (1999,2001), which shall be applied to study 
the planets migration inside a time-dependent accretion disc. Since the model has already been described in the two quoted 
papers, the reader is referred to those for further details. 

We consider a thin planetesimal disc around a star of mass M, — IMq and suppose that a single planet moves in the 
disc under the influence of the gravitational force of the star. The equation of motion of the planet can be written as: 



(Melita & Woolfson 1996), where the term Fq represents the force per unit mass from the Sun, while R is the dissipative 
force (the dynamical friction term-see Melita & Woolfson 1996). In order to take into account dynamical friction, we need a 
suitable formula for a disc-like structure such as the protoplanetary disc. 

We assume that the matter-distribution is disc-shaped and that it has a velocity distribution described by: 



(Hornung & al. 1985, Stewart & Wetherill 1988) where and cr|| are the velocity and the velocity dispersion in the direction 
parallel to the plane while v± and a± are those in the perpendicular direction. We suppose that (7|| and a± are constants and 
that their ratio is simply taken to be 2:1 ( a^^=2a±). 

Then according to Chandrasekhar (1968) and Binney (1977) we may write the force components as: 



1= = Fq + R 



(19) 




(20) 



-B|l«i|| 2\/27rnG^ log Amim2 (mi + m2) 




(21) 



B±vi^ 2V2TvnG^ log Amim2 {mi + 1712) 




(22) 
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where 



Bn = / exp 







^111 1 







Bx ~ / cxp 



and 



2crn + g 2cr2 1 - e2 + g 



1 



[(l + g)2(l-e2 + g)i 



/21 



(23) 



2(jjl + q 2crr^ 1 - e2 + g 



[(l + g)(l-e2 + Q) 



(24) 



3/21 



e = (1 - crl/o-f )"■" (25) 

while n is the average spatial density, mi is the mass of the test particle, m2 is the mass of a field one, and log A is the 
Coulomb logarithm. The frictional drag on the test particles may be written as: 

F = — fc||?7i||e|| — fcj_ui_Lej_ (26) 
where ey and ex are two unit vectors parallel and perpendicular to the disc plane. 

Since damping of eccentricity and inclination is more rapid than radial migration (Ida 1990; Ida & Makino 1992; Del Popolo 
et al. 1999), we deal only with radial migration and we assume that the planet has negligible inclination and eccentricity, 
ip ~ Gp ~ and that the initial heliocentric distance of the planet is 5.2 AU. For the objects lying in the plane, the dynamical 
drag is directed in the direction opposite to the motion of the particle and is given by (Del Popolo et al. 1999): 

F~-fci|i;,|ej| (27) 

In the simulation, we assume that the planetesimals all have equal masses, m, and that m << M, M being the planet 
mass. This assumption does not affect the results, since dynamical friction does not depend on the individual masses of these 
particles but on their overall density. If the planetesimals attain dynamical equilibrium, their equilibrium velocity dispersion, 
(Jm, would be comparable to the surface escape velocity of the dominant bodies (Safronov 1969), while if we consider a two- 
component system, consisting of one protoplanet and many equal-mass planetesimals the velocity dispersion of planetesimals 
in the neighborhood of the protoplanet depends on the mass of the protoplanet. Since the eccentricity is given by; 



en 



20(2m/3MQ)^''^ M < lO^^g 

(28) 

6(Af/3MQ)^''^ M > lO^^g 



(Ida & Makino 1993) where m is the mass of the planetesimals, then the dispersion velocity in the disc is characterized by 
two regimes, being it connected to the eccentricity by the equation: 

(j,n ~ (e^ -f v)^/^Uc (29) 
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where im is the inchnation of planetesimals and «c is the Keplerian circular velocity. In the neighborhood of the protoplanet 
there is a region having a larger dispersion velocity. The width of this "heated" region is roughly given by 4[(4/3)(eJn + im)'J^ + 
12hlia^]^''^ (Ida & Makino 1993) where a is the semi-major axis and hm = ( sm!" )^^"^ ^® "^^^^ radius of the protoplanet. 
The increase in velocity dispersion of planetesimals around the protoplanet decreases the dynamical friction force (see Eq. 
p^ ) and consequently increases the migration time-scale. 

In our model the gas is almost totally dissipated when the planet begins to migrate. We assume that the spatial structure and 
the time evolution of the surface density, Es(_R, t), in planetesimals is described by the disc model introduced in the previous 
section with the further assumption that the initial radial distribution of solids is Tis{t = 0, r) = E(t = 0, r)5, where 5 = 10~^ 
for ice (gas-giant region) and 5 = 6 x 10""^ for silicates (terrestrial planet region), to account for cosmic abundance. 
The total mass in the planetesimal disc, for a fixed viscosity regime, is then: 

R 

Md{t) = / E ■ 2-Krdr 





- «,.(a™ w 

We integrated the equations of motion using the Bulirsch-Stoer method. 

Before going on it is important to discuss a clue assumption of the paper, above introduced, namely that the surface 
density in planetesimals is proportional to that of gas and that the spatial structure and the time evolution of the surface 
density, Es(-R,t), is such that Es(-R, t) oc T,{R,t). 

We know that the evolution of the surface density of gas, Eq. (^), is described by a diffusive-type equation, while that of 
solid particles is an advection-diffusion equation; 



_ 3d_ 
dt r dr 



1 d 

+ -7r 

r or 



Q,kts 

where i^s = the Schmidt number, Sc, is given by: 



(31) 



Sc = (l-^nkts)yi + ^ (32) 

where v is the relative velocity between a particle and the gas, K the turbulent velocity, f2k is the Keplerian angular velocity, 
ts the so called stopping time. If ilkts — > 0, the stopping time is small in comparison with the period of orbital revolution 
and particles are strongly coupled to the gas (particles of size < 1 mm). If fJkis oo, the stopping time is very long in 
comparison with the period of orbital revolution and particles are decoupled from the gas (particles of size > 10'' cm). Then 
the mass distribution of planetesimals emerging from a turbulent disc does not necessarily reflect that of gas. In general, the 
mass distribution of solids evolves due to gas-solid-coupling, coagulation, sedimentation, and evaporation/condensation. After 
© 2000 RAS, MNRAS 000, §-?? 



14 A.Del Popolo et al. 

some preliminary studies of Cassen (1996) and Schmitt et al. (1997), Stepinski & Valageas (1996, 1997) developed a method 
that, using a series of simplifying assumptions, is able to simultaneously follow the evolution of gas and solid particles due 
to gas-solid-coupling, coagulation, sedimentation, and evaporation/condensation for up to lO^yr. The model is based on the 
premise that the transformation of solids from dust to planetesimals occurs fundamentally through the process of hierarchical 
coagulation (other possibility for this transformation are reported in Goldreich & Ward 1973; Barge & Sommeria 1995; Tanga 
et al. 1996). The model is a hypothetical evolutionary scenario that can at least be used to illustrate the difference in time 
evolution between gas and solids (I want to recall that a further element of imprecision is the dust opacitiy which may 
be uncertain in the mm range by a factor of 4-5 (Pollack et al. 1994)). Stepinski's model, goes only to times of 10^ years 
fundamentally because at such long times solids are mostly in large bodies (planetesimals) and Stepinski's method, which 
neglect gravitational interactions between solids, is not longer a reasonable approximation. 

The comparison of the evolution of the gas disc and that of solids is plotted in Fig. 4-5. Fig. 4 shows the evolution of the 
mass of the nebula for a disc of = O.245M0, a = 0.01. The solid line represents the surface density of the gas xlO~^ 
(obtained with our model), while the dashed line the evolution of long lived particles (characterized by evolutionary timescales 
comparable or longer than those of gas) of lO^cm, calculated by Stepinski & Valageas (1996). 

Fig. 5, shows the evolution of the mass of the nebula for a disc of — O.O23A/0, for a = 0.1 (Fig. 5a), a — 0.01 (Fig. 
5b), a = 0.001 (Fig. 5c), a — 0.0001 (Fig. 5d). The dashed line represents the surface density of the gas xO.Ol (obtained with 
our model), while the dot-dashed line the evolution of solids, calculated by Stepinski & Valageas (1997). 

Summarizing, the previous plots. Fig. 4-5, show that the difference in evolution of gas xlO^^ and that of solids becomes 
smaller for smaller values of a. 

Although, as remarked, the evolution of gas and solids is different, I'll use in the present paper the quoted approximation 
Es(-R,t) oc T,{R,t). There are several reasons for this choice: 

a) In the past, but also in several recent works, the study of origin of the Solar System was based upon a steady-state model, 
namely the minimum-mass model (Weidenschilling 1977; Hayashi 1985). Other papers assumed only spatial, but not temporal, 
changes in the nebula and the origin and evolution of the Solar System was studied by means of steady-state accretion disc 
models (see Morfill & Wood 1989; Stepinski et al. 1993). So, the use, in this paper of a time-dependent accretion disc, is surely 
an improvement on the previous papers. 

b) The assumption that the surface density of the planetesimals disc is 0.01 of the gas disc is used by several authors: for 
example, Terquem et al. (2000) (page 3), write: "the surface mass density of the planetesimal disk was derived from E (the gas 
surface density) by noting that in protostellar disks the gas to dust ratio is about 100". This assumption is not "exceptionally" 

© 2000 RAS, MNRAS 000, 0-?? 



Migration of giant planets in a time- dependent planetesimal accretion disc. 15 



used by the prevous authors, but it is almost the "rule" in hterature (e.g., Lin & Papaloizou (1980, page 47); Murray et al. 
(1998)). The paper of Stepinski & Valageas (1996), which is the first attempt to study, numericaUy, the global evolution 
of solids in a protoplanetary discs, recognizes that the above quoted assumption is notheworthy diffused in literature, and 
nowaday it is always used. 

4 RESULTS 

In the present paper, we are fundamentally interested in studying the planets migration due to interaction with planetesimals 
and for this reason we suppose that the gas is almost dissipated when the planet starts its migration. Clearly the effect of the 
presence of gas should be that of accelerating the lose of angular momentum of the planet and to reduce the migration time. 
Similarly to Del Popolo et al. (2001), we assume that the gas disc has a nominal effective lifetime of 10^ years (Zuckerman 
et al. 1995), compatible with several evidences showing that the disc lifetimes range from lO'' yr to 10^ yr (Strom et al. 
1993; Ruden & Pollack 1991). Usually, this decline of gas mass near stars is more rapid than the decline in the mass of 
orbiting particulate matter (Zuckerman et al. 1995). Then the disc is populated by residual planetesimals for a longer period. 
Summarizing our model starts with a fully formed gaseous giant planet of lA/j at 5.2 AU embedded in disc of planetesimals, 
without gas. The model introduced in the previous sections was integrated for several values of the disc surface density or 
equivalently several disc masses: Md = 0.1, 0.01, 0.005, 0.001, 0.0005, 0.0001 Mq (note that we follow the same notation of 
Del Popolo et al. (2001) and so the value of the disc mass, Md, refers to initial gas disc). 

We recall that estimates of the "minimum mass" disc necessary to form our planetary system range from ~ 0.001 to ~ O.IMq 
(Weidenschilling 1977; Boss 1996). Due to the loss of planetesimals ejected through orbital encounters with giant planets (see 
Murray et al. 1998), estimates on the low end of this range may be insufficient to form our solar system; consequently, the 
minimum mass for our solar system may be probably ~ O.OOM© (Boss 1996). The results of the disc model introduced in 
section (2) are plotted in Fig. 1-3. 

Fig. 1 shows the detailed evolution of the surface density and we compare the surface density obtained with the model 
of section (2), solid line, with that of Stepinski (1998), dashed line, for a disc mass of O.IMq, a = 0.01, Tcrit = 1-78 and the 
angular momentum given by: 



R 




(33) 







The angular momentum can be directly obtained using equation (ph, and for a fixed viscosity regime is given by: 



© 2000 RAS, MNRAS 000, §-?? 



16 A.Del Popolo et al. 

'2K{q + l)Bi 



, ^ , -KriJ^oVGAfFo (34) 
2q- ■p + 2 

where 

+ 1 15g + 5-2p 
is the beta function defined as 

which for the disc parameter previously defined is 4 x 10^'^gcm^/s. The plots, from top to bottom, represent the surface density 
at times t — 10^, 10*', lO^yr. The surface density plotted was obtained by imposing, similar to Stepinski (1998), the continuity 
of E at boundaries between different viscosity regimes. The size of the viscosity regimes and the matching of E(r, t) at the r^j 
was obtained as described in section (2). As shown in the figure, our model for E is in good agreement with Stepinski's (1998) 
model. As shown, there are already five different regimes present at t = lO^yr. In general the angular momentum Jd and the 
disc mass Md determine the overall evolution of the nebula and a sets the time scales: lower values of a give rise to longer 
time scales. The number of epochs in the nebula evolution are strictly connected to the initial conditions. Aid and Jd. In fact, 
fixing the initial mass, the smaller the angular momentum, the more concentrated towards the star is the mass distribution. 
Usually there are two or more regimes and consequently different parts of it have a behavior set by the relative viscosity 
regime. Differently from other models, like Dubruelle (1993) and Cassen (1994,1996), the changes in E are not given by a 
simple decay formula (see Stepinski 1998, equation (1)). This is due to the presence of different epochs and since the viscosity 
is described but more than one power-law. Fig. 1 shows that at lO^yr the disc is in the main viscous phase. The accretion 
of mass towards the central star gives rise to a later evolution characterized by a steady decrease in the surface density. The 
evolution proceeds self-similarly. In Fig. 2, we also compared the result of our model, solid line, with the numerical simulations 
of Reyes- Ruiz & Stepinski (1995), dashed line, and E obtained using our model but only the viscosity regime SIGOR, dotted 
line. As shown there is a good agreement between our analytic model and their numerical simulations. The main differences 
between the analytical model and numerical simulations is larger in the innermost regions of the nebula. This is due to the 
fact that in numerical simulations Reyes- Ruiz & Stepinski (1995) has used additional, high-temperature opacity regimes, 
neglected in the analytical model. It is also interesting to note that the analytical prediction for E obtained using only the 
SIGOR regime is also in good agreement with the other two curves at least in the region 0.3-5.2 AU, which is of particular 
interest for the migration of giant planets studied in this paper. Fig. 3 plots the radial distribution of temperature, T{r,t), at 
selected times (lO^yr, lO^yr, lO^yr) for a disc of Aid ~ O.IMq and a — 0.01. The solid line represents T{r,t) obtained with 
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Stepinski's (1998) model, the dashed line the prediction for temperature of the model of this paper and the dotted line the 
approximation given by using only the SIGOR regime. As in the case of E(r, t), the temperature obtained in this paper is in 
good agreement with Stepinski's (1998) model and Reyes-Ruiz & Stepinski (1995). 

In the successive figures (6, 7, 8), we plot the evolution of semi-major axis of the planet. In Fig. 6, we plot the evolution 
of a 1 Mj planet in a disc with planetesimals whose surface density is 1% of the gas (Lin & Papaloizou 1980; Stepinski 
& Valageas 1996) and having a = 0.001. The disc model used is that introduced in section 2. As previously reported, the 
simulation is started with the planet at 5.2 AU and ip ~ ep ~ 0. The values of the initial disc mass are Md: 0.1, 0.01, 0.005, 
0.001, 0.0005, 0.0001 Mq, for the solid line, dotted line, short-dashed line, long-dashed line, dot-short dashed line, dot-long 
dashed line, respectively. As it is evident, a disc of larger mass produces a more rapid migration of the planet. Fig. 6 shows 
that the planet embedded in a disc having Afo = 0.1, 0.01, 0.005, 0.001, 0.0005, 0.0001 Mq migrates to 0.05 AU in 9.5 x lO^yr, 
2.7 X lO^yr, 3.7 x 10*yr, 8 x lO^yr, 1.2 x 10'"*yr, 2.5 x lO^yr, respectively. It is important to stress that even if the planet reaches 
distances of 0.05AU in less than 4.5 x lO^yr, the planet must halt at several i?, from the star surface (i?* is the stellar radius). 
In fact solid bodies cannot condense at distances < 7-R, , and planetesimals cannot survive for a long time at distances < 2_R* . 
When the planet arrives at this distance the dynamical friction force switches ofi' and its migration stops. 
This means that the minimum value of the semi-major axis that a planet can reach is ~ 0.03 AU. 

In order to study the effect of viscosity on migration, we performed two other simulations with a = 0.01 and a — 0.1, and 
plotted in Fig. 7 and Fig. 8, respectively. In Fig. 7, we plot the case a — 0.01. As for Fig. 6 the values of the initial disc 
mass are Mo'- 0.1, 0.01, 0.005, 0.001, 0.0005, 0.0001 Mq, for the solid line, dotted line, short-dashed line, long-dashed line, 
dot-short dashed fine, dot-long dashed line, respectively. The plot shows that for A/d ~ O.IA/0 the planet moves to 0.05 AU 
in ~ 8.8 X 10* yr. If Md = 0.01, O.OO5Af0 , the time needed to reach 0.05 AU is given by ~ 2.5 x lO'' yr, ~ 3.2 x lO'' yr, 
respectively. If the disc mass is A/d = 0.002Mq, the time needed for the planet to reach the quoted position is larger than the 
age of the stellar system. Similar to the case for a = 0.001, even in the cases Mo ~ 0.1, 0.01, O.OO5Af0, in which the planet 
reaches distances of 0.05AU in less than 4.5 x 10® yr, the planet must halt at several R, from the star surface {R, is the stellar 
radius). In the case Mo = O.OOIA^q, the planet stops its migration at a ~ 0.3 AU, while if Mo = 0.0005, 0.0001, Af0, we 
have a ~ 0.7, 1.86 AU. Fig. 8 is the same of the previous two figures except that a = 0.1 in this case. In particular, we find 
that the planet embedded in a disc having Mo = 0.1,0.01, 0.005, 0.001, 0.0005, 0.0001 A/q migrates to 0.53, 2.12, 2.57, 3.43, 
3.73, 4.25 AU, respectively, in 4.5 x lO^yr. We see that the migration time increases going from a = 0.001 to a = 0.1 since 
discs with lower value of a evolve slower than discs with larger values of this parameter, and the high value of the gas and 
particle density is retained for a much longer time (Stepinski & Valageas 1996). As the solar nebula evolution is faster for 
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larger values of a, larger amounts of mass are lost from the nebula, in a given time, and since a less massive disc implies a 
less rapid migration, the planet shall move less, in the case of larger ct (see Ruden & Pollack 1991, Fig.l). 

The final distribution of planets shows that, in the case of discs having a value of a = 0.001, for disc mass in the range 
O.IMq < M < O.OOOIM0, a Jupiter-like planet in any case migrates to a very small distance from the parent star and that, for 
the reason described above, the migration stops at several _R, . If the value of a = 0.01, planets can migrate both to very short 
distances from the star (for disc masses in the range O.IMq < M < O.OO5M0) or to distances ~ 2AU for M ~ O.OOOIM©. 
Planets cannot migrate to very small distances from the star, if a = 0.1 and O.IA/q < M < O.OOOIA^q. In this last case, the 
planet is located in the range 0.5 < a < 4. Summarizing, according to the final distribution of planets distances, similar to Del 
Popolo et al. (2001), the present model predicts that planets can be present at any distance from their locations of formation 
and very small distances from the parent star. Differently from Del Popolo et al. (2001), now the location of the planets also 
depend on the viscosity through the parameter a. It is evident from the previous results that the viscosity has at least the 
same importance of disc mass for what concerns migration. Another difference that our new disc model has on migration of 
planets is the different time scale and that the evolution of the semi-major axis never goes to a constant value as is did in 
Del Popolo et al. (2001). This last difference is due to the more rapid decrease in surface density and disc mass assumed in 
Del Popolo et al. (2001) (see Fig. 2 of Del Popolo et al. (2001)). One of the aims of this last paper was that of finding the 
qualitative efi'ect of disc evolution and this led us to assume a decay law for the disc mass similar to that for cm size particles, 
given in Stepinski & Valageas (1996). This was only a rough assumption, to have an idea of the effect of mass evolution in 
the disc on the migration of planets. In the present paper, the disc model gives us a better background to understand how 
the spatial structure and the temporal evolution of the nebula influences the migration of planets. 

The configuration of observed planets can be reproduced by means of a combination of a and Md- For example, the 
parameters of 47 UMa b (a = 2.11 AU) can be explained, for example, with a — 0.01 and Aid — 0.0001 or a = 0.1 and 
Md — 0.01, while planets like t Bootis b or 51 Peg b, having very small semi-major axis, can be reproduced with a = 0.001 and 
O.IMq < Md < O.OOOIM0 or a = 0.01 and O.IA/q < Md < O.OOSSAfg. Configurations similar to that of 55 Cue h {a = 0.11 
AU), p CrB b (a = 0.23 AU) are obtained for a = 0.01, Md ~ O.OO15M0 and a = 0.01 A'/d ^ 0.0012A/q, respectively. The 
semi-major axis of planets like 70 Vir b (a = 0.43, ep = 0.4; Afsinip ~ 6.6Mj) and HD 114762 b (a = 0.3; Cp = 0.25; 
Msinip ~ lOAfj) can be explained by radial migration, as shown, while their high value of eccentricities can be explained as 
described in Del Popolo et. al (2001). |^ For what concerns the enhancement of metallicity of stars with short-period planets, 

^ If a planet having mass M > 3Afj moves in a planetesimal disc during interactions, planetesimals scattered from their Hill sphere 
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having high metallicities, [Fe/H] > 0.2 (Gonzales (1997; 1998a,b), there is no difference with Del Popolo et al. (2001); namely 
we expect that the planet delivers into the star, in the case of a disc A/d = O.OIM0, a mass of Mace ~ 4OM0 which means 
[Fe/H] ~ 0.2. We can add that, since the quoted mass is delivered if the planet moves to a distance of 0.05AU, the disc should 
preferably have a = 0.001 and Md = O.OIM0, a = 0.01 and A/d = O.OIA^q ov a = 0.1 and Md = Q.IMq (note that this mass 
should be contained in a disc of several tenths of AU). 

5 CONCLUSIONS 

In this paper, we have studied the effect of dynamical friction on the migration of a giant planet in a planetesimal time- 
dependent accretion disc, having only ~ 1% of its mass in the form of solid particles (Stepinski & Valageas 1996, Murray et. 
al 1998). The paper is an improvement of a previous one (Del Popolo et al. 2001), in which the effect of dynamical friction 
has been studied for a disc having surface density E(r) = E0(lAU/r)^'^^, where Sq, the surface density at 1 AU, is a free 
parameter. To start with, we found a self-similar solution to the diffusion equation: the solution obtained is very similar to 
that obtained by Stepinski (1998) in his theoretical model. Our theoretical model for the disc is in good agreement with 
numerical calculations of Reyes- Ruiz & Stepinski (1995) and with the theoretical model of Stepinski (1998). The disc model 
was used to calculate the planet migration. We found in the case of discs having a value of a = 0.001, and for disc mass in 
the range 0.1 A^q < M < O.OOO1A^0, that a Jupiter-like planet in any case migrates to a very small distance from the parent 
star and that, for the reason described above, the migration stops at several i?*. If the value of a = 0.01, planets can migrate 
both to very short distances from the star (for disc masses in the range O.IA^q < Af < O.OO5A'/0) or to distances ~ 2AU for 
M ~ O.OOOlAf0. Planets cannot migrate to very small distances if a = 0.1 and O.lAf0 < M < O.OOOlAf©. In this last case, 
the planet is located in the range 0.5 < a < 4. The location of planets depend on two parameters: the viscosity, through the 
parameter a, and the disc mass. Different configurations of observed planets (e.g., 47 UMa b; r Bootis b; 51 Peg b; 55 Cue b; 
p CrB b; 47 UMa b) can be explained by combinations of a, and the disc mass. The configuration of large eccentricity planets 
like 70 Vir b (a = 0.43, Cp = 0.4; Msinip ~ 6.6Mj) and HD 114762 b (a = 0.3; Cp = 0.25; M sinip ~ 10A//j) can be explained 
by radial migration, while their high value of eccentricities can be explained observing that, in the case M > 3Mj, eccentricity 
the planets ep tends to increase if planetesimals scattered from their Hill sphere can be ejected with j-^j < 1. Similar to Del 
Popolo et al. (2001), metallicity enhancement observed in several stars having extrasolar planets can also be explained by 

can be ejected with |^^| < 1, (where AE and AL are respectively the energy and angular momentum removed from a planet by the 
ejection of a planetesimal), and the eccentricity Ep tends to increase. 
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means of scattering of planetesimals onto the parent star, after the planet reached its final configuration. However, the disc 
should preferably have a = 0.001 and Ma = O.OIMq, a = 0.01 and Md = O.OIMq or a = 0.1 and A/d = O.IMq. 
Interesting points to address in a future study are: 

a) migration in presence of at least two planets; 

b) effect of gas and dynamical friction on the migration; 

c) effect of accretion and change in the mass of the planets. 

Melita & Woolfson (1996) have performed numerical integrations of a three-body problem to study the combined effect of 
resonances and dissipative forces on the production of stable configurations. It should be interesting to improve the model, 
by means of an appropriate disc model, letting mass of the planets change by accretion, and using the formulas for dynamical 
friction force used in this paper, instead of the form they took from the stellar systems evolution theory, and study the effect 
of different choices for the initial masses of the planets. 
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24 A.Del Popolo et al. 

Figure 1. Evolution of tlie surface density. The solid line represents the disc model of this paper (see section (2)) while the dashed 
lino represents Stepinski's (1998) result. In this plot, the disc mass is O.IMq, a = 0.01, t„\^ = 1.78 and the angular momentum is 
4 X lO^^gcm^/s. The plots, from top to bottom, represent the surface density at times t = 10^, 10®, lO^yr. 
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Figure 2. Comparison of tlie disc model of tliis paper, solid line, with same parameters as Fig. 1, with the numerical simulations of 
Reyes-Ruiz & Stepinski (1995), dashed line, and E obtained using our model but only the viscosity regime SIGOR, dotted line. As in 
Fig. 1, the plots, from top to bottom, represent the surface density at times t = 10^, 10®, lO'^yr. 
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Figure 3. Radial distribution of temperature, T{r,t), at selected times (lO^yr, lO^yr, lO'^yr) (from top to bottom) for a disc of 
J\/j = O.IMq, a = 0.01, Tcrit = 1-78 and J^j = 4 X lO^^gcm^/s. The solid line represents T{r,t) obtained with Stepinski's (1998) model, 
the dashed line the prediction for temperature of the model of this paper and the dotted line the approximation given by using only the 
SIGOR regime. 
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Figure 4. Evolution of tlie mass of the nebula for a disc of = 0.245Mq, a = 0.01. The solid line represents the surface density of 
the gas xO.Ol, while the dashed line the evolution of long lived particles (see Stepinski & Valageas 1996) of 10*cm. 
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Figure 5. Evolution of the mass of the nebula for a disc of A/j = O.O23M0, for a = 0.1 (a), a = 0.01 (b), a = 0.001 (c), a = 0.0001 
(d). The dashed line represents the surface density of the gas xO.Ol, while the dot-dashed line the evolution of solids. 
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Figure 6. Evolution of semi-major axis of a 1 Mj planet in a disc with planetesimals whose surface density is 1% of the gas and 
having a = 0.001. The values of the initial disc mass are Mq: 0.1, 0.01, 0.005, 0.001, 0.0005, 0.0001 Mq, for the solid line, dotted line, 
short-dashed line, long-dashed line, dot-short dashed line, dot-long dashed line, respectively. 
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Figure 7. Same as Fig. 4 but now a = 0.01. 
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